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Abstract We use probability density functions (pdfs) of sums of orbit coordinates, over time intervals of 
the order of one Hubble time, to distinguish weakly from strongly chaotic orbits in a barred galaxy model. 
We find that, in the weakly chaotic case, quasi-stationary states arise, whose pdfs are well approximated 
^ — ■ by g-Gaussian functions (with 1 < q < 3), while strong chaos is identified by pdfs which quickly tend to 

Gaussians (q = 1). Typical examples of weakly chaotic orbits are those that "stick" to islands of ordered 
motion. Their presence in rotating galaxy models has been investigated thoroughly in recent years due of 
their ability to support galaxy structures for relatively long time scales. In this paper, we demonstrate, on 
specific orbits of 2 and 3 degree of freedom barred galaxy models, that the proposed statistical approach 
can distinguish weakly from strongly chaotic motion accurately and efficiently, especially in cases where 
pH , Lyapunov exponents and other local dynamic indicators appear to be inconclusive. 
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■ 1 Introduction 

o i 

In galactic dynamics, the identification of orbits as ordered (i.e. periodic or quasiperiodic) or chaotic 
and their respective role in the structure of galaxies has been studied extensively for many years (see 
|12] for a recent review). Distinguishing between quasiperiodic and chaotic motion, in general, poses no 
major difficulty, as there are many numerical approaches presently available by which this distinction 
can be easily accomplished. What is considerably more complicated is the identification of weakly chaotic 
orbits, particularly when integration times are not long enough. Weak chaos is manifested, for example, 
in trajectories that "stick" for long times on the boundaries of "islands" of ordered motion, where the 
maximal Lyapunov exponent (MLE) is positive but "small" , before eventually occupying a wider domain 

■ of strong chaos with much larger MLEs. 
Weakly chaotic orbits spend a significant fraction of their time in confined regimes and do not fill up 

phase space as "homogeneously" as strongly chaotic trajectories, which explore faster and more uniformly 
their available regions. Several researchers have demonstrated the importance of weakly chaotic phenomena 
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in supporting barred or spiral galaxy features (sec e.g. 3,4 518,21 30 , 34 , 40 , 41 ). There are also several 
results in the recent literature showing that strong local instability does not necessarily imply widespread 
diffusion in phase space [10II11II17] . In [131114] "stickiness" was studied thoroughly in 2 degrees of free- 
dom (dof) and was categorized in two main types: (a) "Stickiness" around an island of stability and (b) 
"stickiness" in chaos, along the invariant manifolds of unstable periodic orbits. The role of "sticky" chaotic 
orbits and the diffusive behavior in the neighborhood of invariant tori surrounding periodic solutions of 
the Hamiltonian was studied in [23]. The challenging problem, however, is to be able to detect weak chaos 
(as manifested e.g. in "stickiness" phenomena), for short integration times (not exceeding the age of the 
universe), before the true character of the dynamics becomes apparent. 

In the present paper, we propose a novel approach for tackling this problem based on the statistics 
of the orbits of certain barred galaxy models. In particular, we consider probability density functions 
(pdfs) of sums of their variables, in the spirit of the Central Limit Theorem If these pdfs quickly 
tend to Gaussians, we classify the orbits as strongly chaotic and suggest that their motion is characterized 
Boltzmann-Gibbs statistical mechanics. If, on the other hand, the pdfs are well approximated by g-Gaussian 
functions, we characterize the orbits as weakly chaotic, implying that they form quasi-stationary states 
that obey the principles of nonextensive statistical mechanics 1 , 43 

Pdfs of chaotic trajectories of dynamical systems have long been studied by many authors, aiming to 
understand the transition from deterministic to stochastic dynamics [2, 15, 22 ,35 ,42, 45]. The main ques- 
tion here concerns the existence of an appropriate invariant probability measure, characterizing chaotic 
phase space regions where solutions generically exhibit exponential divergence of nearby trajectories. If 
this invariant measure turns out to be a continuous and sufficiently smooth function of the phase space 
coordinates, one can invoke the Boltzmann-Gibbs microcanonical ensemble and attempt to evaluate all 
relevant quantities of equilibrium Statistical Mechanics, like partition function, free energy, entropy, etc. 
of the system. 

In such cases, viewing the values of one (or a linear combination) of components of a chaotic solution 
at discrete times t n ,n = 1, ...,J\f as realizations of M independent and identically distributed Gaussian 
random variables X n and calculating the distribution of their sums, one expects to find a Gaussian distri- 
bution with the same mean and variance as the X n 's (according to the Central Limit Theorem). This is 
indeed what happens for many chaotic dynamical systems studied to date which are ergodic, i.e. almost 
all their orbits (except for a set of measure zero) pass arbitrarily close to any point of the constant energy 
manifold, after sufficiently long enough times. In these cases, at least one of the Lyapunov exponents [BE 
15,50 is positive, stable periodic orbits are absent and the constant energy manifold is covered uniformly 
by chaotic orbits, for all but a (Lebesgue) measure zero set of initial conditions. 

In the galactic models treated in this paper, we focus on weakly chaotic regimes and demonstrate by 
means of numerical experiments that pdfs of sums of their variables do not rapidly converge to a Gaussian, 
but are well approximated, for long integration times, by the so-called g-Gaussian distribution [43] 



where the index q satisfies 1 < q < 3, /3 is an arbitrary parameter and a is a normalization constant. 
At longer times, of course, chaotic orbits eventually "seep out" from small regions to larger chaotic seas, 
where obstruction by islands and cantori is less dominant and the dynamics is more uniformly ergodic. This 
transition is signaled by the g-index of the distribution (|14[) decreasing towards q = 1, which represents 
the limit at which the pdfs become Gaussian. 

In what follows, we begin in Section [5] by presenting our galaxy barred galaxy model with all its 
components and parameters. In Section [3] we describe the approach of the g-Gaussian distributions and 
discuss the details of their computation. Sections [4] [5] are devoted to the application of these pdfs and 
the demonstration of their effectiveness in distinguishing between strong and weak chaos in 2 and 3 dof 
cases of the galactic potential, using also power spectra to support our findings. Finally, we conclude with 
a summary and discussion of our results in Section [6] 

2 The barred galaxy model 

Let us consider the motion of test particles in a rotating "mean field" model of a barred galaxy potential 
V(x, y, z) governed by the 3 dof Hamiltonian 



P(s) = ae X p q (-Ps 2 ) 



a 1 - (1 - q)ps 2 



l-q 



H = 




(1) 
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The bar rotates around its z-axis (the short axis), while the ^-direction rotates along the major axis and 
the y along the intermediate axis of the bar. p x , Py and p z are the canonically conjugate momenta, fif, 
represents the pattern speed of the bar and Ej is the total constant energy of the orbit in the rotating 
frame of reference (namely the Jacobi integral). The corresponding equations of motion are 

x = p x + O b y, y = p y -Q b x, z = Pz, (2) 

dV _ dV _ dV 

Px = --%- + i*bPv} Py = ~~^ "bPx, P z = ~-^~- 

ox ay oz 

The potential consists of three components V = Vrj + V$ + Vg: 

1. The disc, which is represented by a Miyamoto-Nagai disc |29] : 

Vp = - , GM ° (3) 



; 2 + J/ 2 + (A + Vz 2 +B 2 ) 2 

where Mr) is the total mass of the disc, A and B axe its horizontal and vertical scale-lengths, and G is 
the gravitational constant. 

2. The bulge, which is modeled by a Plummer sphere [37J whose potential is: 

v 3 =- /2 G 2 Ms 2 , ■ w 

\/x l + y l + z l + ef 

e s is the scale-length of the bulge and Mg is its total mass. 

3. The triaxial Ferrers bar [16] . the density p(x) of which is 

J Pc (l-m 2 ) 2 , m< 1 

PW = \ n ... (5) 

I 0, m > 1 

where p c = G abc IS ^ ne cen tral density, Mg is the total mass of the bar and 

2 2 2 

m = % + \rr + 2=., a>6>c>0, (6) 
a z o z c z 

with a, fo and c being the semi-axes. The corresponding potential is: 

pOO j 

Pc I alt M 2/ \\n+l 



V B = -*Ga b c^- 1 ^ m (l-m<(u)T^, (7) 
where 

2 2 2 

rn(u) = h -pf h -n , (8) 

w a 2 + u b 2 + u c 2 + u w 

A 2 (u) = (a 2 + it) {b 2 +u)(c 2 + it), (9) 

n is a positive integer (with n = 2 for our model) and A is the unique positive solution of 

m 2 (A) = 1, (10) 

outside of the bar (m > 1) and A = inside the bar. The corresponding forces are given analytically 
in [SH]. 

Throughout the paper, we use the following parameters: G=l, J?b=0.054 (54 km ■ sec -1 • fcpc -1 ), a=6, 
b =1.5, c=0.6, A=3, S=l, e s =0.4, M B =0.1, M s =0.08, M D =Q.82, both for the 2 and 3 dof versions. The 
units used are: 1 kpc (length), 1000 km ■ sec -1 (velocity), 1 Myr (time), 2 x 10 11 A/q (mass). The total 
mass G(Mg + Mjj + Mg) is set equal to 1 and the corotation radius is equal to R c = 6.13. 

The above model was used extensively in orbital studies [31 132, 33. 47) 48]. in which the stability of 
its main periodic orbits was thoroughly investigated, as several of its parameters were varied. In [25,26 
I27| . a dynamical study of regular, chaotic and weakly chaotic motion was presented (both in phase and 
configuration space) for different distributions of initial conditions and values of the parameters of the bar. 
Moreover, the SALI/GALI method 46,49 and frequency spectrum analysis were applied to discriminate 
between strongly and weakly chaotic orbits in 24,28 . The latter are of great importance from an astro- 
nomical point of view, since they can last for long time scales and behave quite "regularly" before revealing 
their inherent chaotic nature. 

In our study, the maximal integration of the orbits is typically set to tt — 10000 Myr (10 billion years), 
which corresponds to a time less than (but of the order of) one Hubble time. Nevertheless, we often extend 
our simulations to longer times to study the dynamics of our solutions in greater detail. 



i 



3 Statistical distributions of chaotic quasi-stationary states and their computation 

The approach we shall adopt is in the spirit of the Central Limit Theorem: Using the solutions of the 
equations of motion (J2j we construct probability density functions (pdfs) of suitably rescaled sums of an 
observable function % = f?(ti), which depends linearly on the position coordinates of the solution and is 
calculated at the ends of M successive time intervals At = t, — tj_i, i = 1, . . . , M. If the motion is 
chaotic and At is relatively "large", the r/(ti) can be assumed to be independent and identically distributed 
random variables (in the large M limit). Thus, evaluating their sum 

M 



s# = £* w (id 



i=l 

for j = 1, . . . , Ni c initial conditions, we can study the statistics of these variables centered about their 
mean value (S«') and rescaled by their standard deviation a^j as 

w - A- - (sg>>) = _L ( f; tf> - ^ ^ -«) 



where 



£ 

3=1 

Plotting the normalized histogram of the probabilities P(s^)) as a function of s^, we approximate the 
resulting pdfs by g-Gaussian functions of the form 



P(s) = aexp 9 (-^s 2 ) 



(1 - q)Ps 2 



(14) 



where q is the so-called entropic index, f3 is a free parameter and a a normalization constant [43]. If the 
motion is strongly chaotic, the pdfs are expected to converge to the well-known Gaussian, which represents 
the limit of (|14[) as q — > 1. However, in regimes of weak chaos it has been found in many studies that 
the orbits form quasi-stationary states, which can be very long-lived and whose pdfs are well-fitted by 
g-Gaussian functions with q well above unity fni6ll7Tf43] . It can be shown that (|14[) is normalized by setting 



r 



3-q 

2(9-1) 



(15) 

(r is the Euler Gamma function), implying that the allowed values of q are 1 < q < 3. The index q is 
connected with the Tsallis entropy [43] 

1-T W o q W 

Sq = k ^ With J2 Pi = 1 (16) 

q 1=1 

where i = 1, . . . ,W counts the microstates of the system, each occurring with a probability and k is the 
so-called Boltzmann universal constant. Just as the Gaussian distribution represents an extremal of the 
Boltzmann-Gibbs entropy Sbg = — k Ej=i Pi m Pii so i s the g-Gaussian (|14p derived by optimizing 



the Tsallis entropy (|16[1 under appropriate constraints. 

Systems characterized by the Tsallis entropy occur very frequently in a variety of applications and obey 
statistics different than Boltzmann-Gibbs, in the sense that their entropy is nonadditive and generally 
nonextensive [43|I44| . Examples falling in this class are physical systems governed by long range forces, 
like self-gravitating systems of finitely many mass points, interacting black holes and ferromagnetic spin 
models, in which correlations decay by power laws rather than exponentially [43] . 

Let us now describe the numerical approach we follow to calculate the above pdfs. First we select rj(t) 
as one (or a linear combination) of the position coordinates x(t),y(t), z(t) of an orbit which we wish to 
identify as strongly or weakly chaotic. Choosing then an integration interval < t < tf sufficiently long to 
obtain reliable statistics, we divide tf into N^ c equally spaced, consecutive time windows, which are long 
enough to contain a significant part of the orbit. Next, we subdivide each such window into a number M 
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of equally spaced subintervals and calculate the sums JTT} at the ends of these subintervals as described 
above. 

In this way, we treat the point at the beginning of every time window as a new initial condition and 
repeat this process Ni c times to obtain as many sums as possible for the time interval tf. Consequently, at 

the end of the integration, we evaluate the Ni c rescaled quantities s^j ((121) and plot the histogram P(s^) 
of their distribution. If we are in a regime where these distributions are well-fitted by g-Gaussians for fairly 
long times, the orbit can be classified as weakly chaotic. If, on the other hand, the orbit is strongly chaotic, 
its pdf is seen to approach a Gaussian, already for relatively short time intervals not exceeding, in our case, 
the age of the universe. 



4 Statistical distributions of the 2 dof model 

Let us start by analyzing the 2 dof version of our galactic model {TJ, setting z(0) — p z (0) = initially, 
whence z(t) = Pz(t) = for all t. In particular, we are interested in identifying the dynamical behavior of 
the following two orbits: 

(i) The 2-Dimensional Strongly Chaotic (2DSC) orbit, with initial condition 
(x, y,p x ,Py) = (0, -0.625, -0.314512, -0.24) and Ej = -0.3, 

and 

(ii) The 2-Dimensional Weakly Chaotic (2DWC) orbit, with initial condition 
(x, y,p x ,Py) = (0, -0.625, -0.002, -0.24) and Ej = -0.36. 

These orbits for the 2 dof case (and the next two of section 5 for the 3 dof case) are selected from a vast 
number of trajectories (see [24, 27]), as representative examples of weakly and strongly chaotic motion to 
which we can apply our pdf approach and link the results with questions of interest to dynamical astronomy. 
The authors in [24, 27] used dynamical chaos detectors on thousands of orbits, in order to investigate and 
chart the fractions of chaotic and regular regimes of the phase space and then link these fractions to the 
main parameters of the bar component. Since it was evident from these studies that there is a significant 
amount of weakly chaotic orbits that appear "regular" for astronomically bounded time scales, we decided 
to revisit this question in the present paper from a statistical point of view. 

In Fig. [1] the left panels correspond to the integration time tt = 10000 and the right panels to tf = 
100000. Figs. [T] (a) and (b) represent surface of section plots (y,Py) with p x > 0, x = for the orbit 2DSC, 
with Ej — —0.3, showing clear signs of widespread chaos as the orbit "fills out" rather uniformly a wide 
chaotic domain. Figs.[T](c) and (d), on the other hand, depict the numerically computed pdfs (solid curve) 
for the observable r\ — x + y and the Gaussian distributions (dashed curve) that best fit the data. As 
is evident from these results the proximity of the computed pdfs to a Gaussian, already at tt — 10000, 
suggests that the 2DSC orbit can be classified as strongly chaotic. 

Let us now turn to our analysis of orbit 2DWC in Fig. [5] with Ej = —0.36 using again integration 
times tf — 10000 for the left panels and tf = 100000 for the right panels. Note on the surface of section 
panels of Fig. EJa) and (b) that the 2DWC orbit fills a region that is considerably more limited than that of 
the 2DSC orbit (see Figs. H>),(b)), containing several islands to which the orbit apparently sticks for long 
times. Indeed, when we compute the sum pdfs of the variable rj = x + y for this orbit in Figs.[!3c) and(d), 
we find that they are well-approximated by g-Gaussians that are quite far from a Gaussian, already for 
tf — 10000, thus suggesting that 2DWC cannot be regarded as a strongly chaotic orbit. Note, however, 
that "stickiness" is just one aspect of this form of chaos, since the 2DWC orbit spends only part of its 
evolution in the vicinity of a number of major "islands" . 

It is important to note, however, that the distribution functions plotted in Figs.[5Jc) and (d) have q > 3 
and thus are not normalizable and cannot be used to represent g-Gaussians for this orbit. They only serve 
to show that the dynamics is certainly far from Gaussian and hence the 2DWC orbit cannot be classified 
as strongly chaotic. 

In Fig. [3]we plot the time evolution of the MLE of the two orbits of the 2 dof system. We find that the 
MLE of the strongly chaotic chaotic orbit 2DSC (solid curve) converges to the value a ~ 0.03011, while 
the MLE of the weakly chaotic orbit 2DWC (dashed curve) tends to a ~ 0.01489. Although this value is 
quite smaller than the one for the 2DSC orbit, these results demonstrate that merely knowing the values 
of the MLEs is not sufficient for determining whether an orbit belongs to the weakly or strongly chaotic 
family. 




Fig. 1 Both rows correspond to tt = 10000 (left panels) and tt = 100000 (right panels), (a) and (b): surface of 
section plots (y,p y ) with p x > 0, x = for the 2DSC orbit of the 2 dof Hamiltonian system ifjQl, showing evidence of 
widespread chaos (the horizontal scale is in kpc in both top panels), (c) and (d): Linear-log scale plots of the pdf of 
the numerically computed orbit (solid curve) and a Gaussian distribution (dashed curve) fitted to the data for the 
observable r\ = x + y. In (c) we use JVj c = 4000 time windows and M = 50 terms in the computation of the sums and 
the numerical fitting with the ^-Gaussian i l 1 41 ) gives q 1.095 with x 2 ~ 0.0003. (d) corresponds to N lc = 20000 
and M = 10 and gives q 0.97 with x 2 ~ 0.00033. Since in both panels (c) and (d) the pdfs are well approximated 
by Gaussians, we conclude that 2DC is a strongly chaotic orbit 



5 Statistical distributions of the 3 dof system 

Let us now focus on the 3 dof system ([T]) and carry out a similar study as in the previous section for the 
orbits: 

(iii) The 3-Dimensional Strongly Chaotic (3DSC) orbit, with initial condition 
(x,y,z,p x ,p y ,p z ) = (0.5875,0,1.291670,0,0,0) with Ej = -0.2792149022676664, 
and 

(iv) The 3-Dimensional Weakly Chaotic (3DWC) orbit, with initial condition 
(x,y,z,p x ,p y ,p z ) = (2.35,0,0.08883,0,0.133330,0) and Ej = -0.2852654501087482 

of the 3 dof Hamiltonian system {1} . 

Figs. HJa) and (b) show projections of the orbit 3DSC in the (x,y) and (x,z) planes respectively for 
the integration time interval tf = 10000. Evidently, no sign of any regularity is observed here and the 
plots indicate that the orbit might be strongly chaotic. Indeed, as we see in Figs. [If c) and (d), both for 
tf — 10000, pdf plots of the variables r\ = x + y in (c) and r\ = z in (d) provide evidence of strong chaos, as 
both distributions are well-approximated by Gaussians, whose q is close to unity. Panels (e) and (f) show 
pdf plots for the variables r\ = x + y and z respectively but now for integration time tt = 100000. Here 
also, the numerically fitted curves are close to Gaussians. 

Note, however, that the upper part of the numerical data in Fig.[4ff), is also approximated rather well 
by a g-Gaussian with q ss 1.94, which indicates a peculiar property of the z projection of the dynamics, 
which will be discussed in detail below. Still, based on these findings, we suggest that the 3DSC orbit be 
classified as strongly chaotic. 

Performing now a similar study for the 3DWC solution, we present in Figs. [5fa),(b) two projections 
of the orbit for tf = 10000, on the (x, y) and (a;, z) planes respectively. This orbit is very close to the so- 
called "xl" family of periodic orbits, which have a barred morphological shape and are generally unstable 
(see [17] , Table 2 for more details). A certain regularity does appear here, especially in the [x, y) plane of 
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Fig. 2 Wo use again final integration times tf = 10000 for the left panels and tf = 100000 for the right ones, (a) 
and (b): surface of section plots of {y,p y ) with p x > 0, x = for the orbit 2DWC of the 2 dof galactic system (TjQ| 
(the horizontal scale is in kpc in both top panels), (c) and (d): Linear- log scale plots of numerical (solid curve), 
g-Gaussian (dotted curve) and Gaussian distributions (dashed curve) for the observable rj = x + y, showing evidence 
of weak chaos. In (c) we have used Ni c = 4000 time windows and M = 50 terms in the computation of the 
sums and the numerical fitting with the <j-Gaussian H14I) gives q f» 3.52 with x 2 ~ 0.00074. In (d) we have taken 
N ic = 25000 and M = 50. Here, the numerical fitting gives q w 3.539 with x 2 ~ 0.00057. Of course, the fact that 
the entropic parameter is q > 3 implies that these distribution functions arc not normalizable and hence cannot 
represent q-Gaussians for this orbit. 




t 



Fig. 3 Time evolution of the MLE of the strongly chaotic orbit 2DSC showing a tendency to converge to a ~ 0.03011 
(solid curve) and the MLE of the weakly chaotic 2DWC orbit approaching the value a ~ 0.01489 (dashed curve), 
for the 2 dof Hamiltonian system (I). 

panel (a). Both panels suggest that the motion predominantly occurs in that plane, while frequent chaotic 
"excursions" are also visible in the z direction. In Figs. [SJc),(d) we plot the numerical (solid curve), q- 
Gaussian (dotted curve) and Gaussian (dashed curve) pdfs of the 3DWC orbit. Fig. [5jc) corresponds to 
rj — x + y, for integration time tf = 100000, which still indicates a good approximation by a g-Gaussian 
with q w 2.464, suggesting that 3DWC is weakly chaotic. 

Interestingly enough, the same cannot be said for rj — z in Fig.[S{d), where the data at tt = 10000 "look 
like" a Gaussian but are not statistically reliable to draw any conclusion. To investigate the z dynamics 
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Fig. 4 Projections of the orbit 3DSC of the 3 dof Hamiltonian system for integration time tj = 10000 in (a) the 
(x,y) plane and (b) the (x, z) plane (the axes' scales are in kpc in both top panels). Second row: Plots in linear-log 
scale of numerical (solid curve), q-Gaussian (dotted curve) and Gaussian (dashed curve) distributions of the same 
orbit, (c) for the variable rj = x + y, with N ic = 4000, M = 50 and (d) r] = z, for N ic = 4000 and M = 50, both for 
tf = 10000. Here, the numerical fitting with I114H gives q as 1.25 with \ 2 ^ 0.00017. In the third row we have used 
tf = 100000 and plotted pdfs for: (e) the variable i) = x + y, N{ c = 10000 and M = 100, fitted by a (/-Gaussian with 
q 0.95 and x 2 ~ 0.00029, and (f) the r) = z variable, with N ic = 10000, M = 100, where the numerical fitting 
parameters are q Rj 1.94 and x 2 ~ 0.00035. 



further, therefore, we plot in Fig. EJe) its pdf for longer integration times (tf = 100000) and discover that 
its central part is well-fitted by a (/-Gaussian, with q w 2.262 with \ 2 ~ 0.000325, suggesting that some 
type of "order" is now present. Note also that, in comparison with the 3DSC orbit (see Fig. Hff)), the 
statistical distribution of the z-coordinate has a higher q value and thus we conclude from all available 
evidence that 3DWC is a weakly chaotic orbit. 

Finally, in Fig. [6l we present the time evolution of the MLE of the strongly chaotic orbit 3DSC, which 
does appear to converge rather rapidly to the value a ~ 0.02600 (solid curve). In the case of the weakly 
chaotic orbit 3DWC, we observe a much slower convergence as the MLE "dives" to small values w 3 x 10~ 3 
(dashed curve) at about tt « 5000, indicating a weaker form of chaos over this interval and only rises to 
higher values corresponding to strong chaos for longer times. In fact, in the time interval < t < 3000, 
the 3DWC orbit does display some quasiperiodic features (see also Figs[5£a),(b)), which will be discussed 
below, when we study amplitude spectra of the 3DWC orbit. 
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Fig. 5 First row: Projections of the orbit 3DWC of the 3 dof Hamiltonian system £Q| for final integration time 
tf = 10000 (a) on the (x,y) plane and (b) on the (x, z) plane (the axes' scales are in kpc in both top panels). Second 
row: Plot in linear- log scales of numerical (solid curve), g-Gaussian (dotted curve) and Gaussian distributions (dashed 
curve) for tj = 100000. (c) corresponds to r\ = x + y, 7V; C = 20000, M = 100 and presents numerical fittings with a 
g-Gaussian (13). yielding q « 2.464 with \ 2 ~ 0.00101. In case (d), however, which corresponds to r] = z, Ni c = 4000 
and M = 50 the data are insufficient to provide reliable estimates. Studying the variable r\ = z for longer times 
tf = 100000, with Ni c = 20000 and M = 100, a satisfactory numerical fitting by a q-Gaussian is obtained giving 
q Ri 2.262 with x 2 ~ 0.00033. 



Intrigued by the above results, we decided to investigate separately the projections of the 3DWC orbit 
on the (x, y) plane and the z axis, for longer integration times. In particular, in Fig. [7]we present the time 
dependence of the entropic index q averaged over M intervals, {q)M, taking M = 50, 100, . . . , 450, 500 and 
N ic = 20000 time windows, computed over intervals [0, T] where T = 10000, 20000, . . . , 990000, 1000000. 
The mean and standard deviation of (q)M was found for all pairs of (A^ c , M) to yield normalized pdfs, i.e. 
1 < <1 < 3 for each time window [0, T]. 

Thus, in Fig. [TJa) we have plotted (<?)m for the observable rj — x + y (see Fig. E{ a )) an d in Fig. [TJb) 
for rj — z (see Fig. Ob)). For r\ = x + y, there is quite a long transient time (nearly 85000 units), for which 
the motion on the (x, y) projection gives no statistical data (i.e. no pair of (iVj c , M) that yields 1 < q < 3), 
since the motion displays strong quasiperiodic features in that time interval. This is also evident by the 
fact that the MLE of the orbit (dashed curve in Fig. [S| decreases on the average to fairly low values over 
this interval. A statistically reliable fit is possible only after about t ~ 90000, yielding values well 



10 




10 1 10 2 10 3 10 4 10 5 



t 

Fig. 6 The time evolution of the MLE of the chaotic orbit 3DSC (it ~ 0.02600: solid black curve) and of the weakly 
chaotic orbit 3DWC of the 3 dof Hamiltonian system JQl (a ~ 0.0236: dashed black curve). 



above 2 with a trend to approach q = 1 at longer times (see Fig. [7|. This is in compatible with our claim 
that 3DWC is a weakly chaotic orbit, for short times. 

In the case of the rj = z observable however, the opposite happens: As we see in Fig. [Tfb) , in the 
beginning, (q) m shows a tendency to diminish towards q = 1 at about t = 30000, although the statistics 
is unsatisfactory, since most of the (Ni c ,M) pairs yield q < 1 or q > 3. Soon after, however, it increases 
significantly until it reaches a maximum q ~ 2.28 at t — 130000 before it starts to decrease again to 
q w 1.16 at t — 1000000. These results are in agreement with our earlier observation that in the z direction 
the dynamics is initially more strongly chaotic than in the x, y projections, but does show a clear tendency 
towards weak chaos at later times. 




10 4 10 5 10 s 10 4 10 5 10 6 

t t 

Fig. 7 The time evolution of the averaged entropic parameter (q)m (see text) for the weakly chaotic orbit 3DWC 
of the 3 dof system Q for (a) r) = X + y, N ic = 20000 and M = 50, 100, . . . , 450, 500. (b): Same as in panel (a) for 
rj = z. The dashed lines correspond to one standard deviation from the average entropic parameter. Note that the 
horizontal axes are in logarithmic scale and the vertical in linear scale. 



This peculiar behavior of the z dynamics is better understood when we compute amplitude spectrum 
of the z component of the motion, using Fast Fourier Transform (FFT) techniques [38]. Thus, in Fig. [8] 
we present such spectra for 4 time windows as follows: In Fig. [8ja) the spectrum is shown for the initial 
time window [0, 10000], where the amplitude of the main frequency u) z « 0.07 is approximately 8, 500. In 
Fig. |5|[b), the same plot is presented for the time window [40000, 50000], with the amplitude of the main 
frequency increasing to a value « 19000 (in accordance with the increase of the (g)jvf parameter found for 
this interval). 

In Fig. [Sfc), the above spectrum is plotted for the time window [90000, 100000], where the amplitude 
of the main frequency now increases to about 27000 (in this time window (q)m nas increased even more, 
reaching a value about 2.1). Finally, in Fig. [BJd) we plot the spectrum for t £ [990000, 1000000] and see 
that as time increases, the energy of the spectrum is restricted within a very limited set of frequencies 
around a central ui z value, which seems to continuously gain more energy and play a dominant role that is 
most likely responsible for the weakly chaotic behavior of the z component of the orbit. 
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Fig. 8 Amplitude spectrum of the coordinate component z for the weakly chaotic orbit 3DWC of the 3 dof system, 
for 4 time windows of length 10000. (a) Shown here is the spectrum for the time window [0, 10000] where the ampli- 
tude of the main frequency u) z Rj 0.07 is approximately 8500. (b) Same as in (a) for the time window [40000, 50000] 
where the power of the main frequency has now increased to a value ~ 19000, a trend that agrees with the increase 
of (q)M over this time interval, (c) Spectrum is plotted for the time window [90000, 100000], where the amplitude 
of the main frequency has further increased to a value of about 27000, while (q)M increases even more to a value of 
about 2. Finally, in panel (d) we plot the same spectrum for t S [990000, 1000000]. 
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Fig. 9 Amplitude spectrum of the coordinate component x + y for the weakly chaotic orbit 3DWC of the 3 dof 
system, for 4 time windows of length At = 10000. (a) The spectrum is shown for the time window [0, 10000] where 
the amplitude of the main frequency ui x +y ~ 0.02 is approximately 150000. (b) Same as (a) for the time window 
[40000, 50000] , where the dominant frequency is now almost absent and the total spectrum has spread to a denser 
set of frequencies. In panel (c), we plot the spectrum for the time window [90000, 100000]. It is in this interval, that 
(<?)m starts gradually to decrease to 1, the limit at which the distribution becomes Gaussian. Finally, in panel (d) 
we plot the spectrum for t £ [990000, 1000000]. 

To provide further support for our claim that 3DWC is weakly chaotic, we performed the same kind 
of analysis for its x + y component in Fig. [9] Fig. [9ja) shows the frequency spectrum for the time window 
[0, 10000], where the amplitude of the main frequency 0J x +y ~ 0.025 is nearly 150000. In Fig. [9|b), the 
spectrum is presented for the intermediate time window [40000, 50000] , over which the data has spread to 
a wider set of frequencies. In Fig. [9jc), the spectrum for the time window [90000, 100000] look similar to 
panel (b). Interestingly, over this time interval, (q) m still increases, reaching its maximum value q ~ 2.28 
at t = 130000. 

Fig. [5£d), on the other hand, for t £ [990000, 1000000] clearly shows that a broad set of low amplitude 
frequencies still carries the energy of the orbit, without any dominant frequency present, of the type found 
in Figs. |[a)-(c). This further demonstrates that the orbit becomes more strongly chaotic only after very 
long times, as shown also by the approach of the (q)M index to q = 1 in Fig. E|a,)- 



6 Summary and conclusions 

We have studied statistical distributions of sums of position coordinates for some typical examples of 
chaotic orbits of a 2 and 3 dof barred galaxy model, aiming to classify them as strongly or weakly chaotic 
over relatively short time intervals (of the order of one Hubble time). This classification is important, when 



12 



one wishes to quantify order and chaos in galaxy models, and is not easy to achieve for short times using 
only local dynamic indicators such as Lyapunov exponents. 

More specifically, we have analyzed the statistics of two orbits of a 2 dof barred galaxy potential, having 
z(t) =0 for all time, and two orbits of the 3dof model, whose z(t) components evolve away from zero. If 
the pdfs produced by the orbits over these intervals are well-approximated by (/-Gaussian functions, with 
q significantly larger than 1 (1 < q < 3), we identify the orbits as weakly chaotic. If, on the other hand, q 
is close to 1, the pdf is well-fitted by a Gaussian and the orbits are said to exhibit strong chaos. 

In summary, we have arrived at the following conclusions: 

1. Weak chaos is a concept that encompasses more that the well-known "stickiness" phenomenon observed 
near the boundaries of major "islands" of ordered motion or along invariant manifolds of unstable 
periodic orbits. 

2. Statistical methods based on q-Gaussian (1 < q < 3) distribution functions can reliably distinguish 
between strongly and weakly chaotic dynamics in barred galaxy models, over relatively short time 
intervals. 

3. In general, weakly chaotic orbits evolve into wider chaotic domains with Gaussian pdfs corresponding 
to strong chaos, but over longer time scales exceeding one Hubble time. This shows that weak chaos 
is an important feature of stellar dynamics that may possibly influence the formation and support of 
certain galaxy structures. 

4. Statistical distributions of different observables do not always show the same dynamical features in 
different projections. For example, the pdfs of the rj = x + y and r\ — z variables of the galaxy model 
studied here, revealed a weaker form of chaotic motion (and for longer time intervals) in the x,y plane 
compared to the z direction. 

5. The averaged entropic index (q)m is a useful parameter for monitoring the nature of the dynamics of 
different observables. In the case of our 3 dof weakly chaotic orbit for example, the surprisingly distinct 
time evolution of (q) m f° r the r\ — x + y observable, compared with T) = z, reflects the sharply different 
dynamics of the orbit in the x, y and z projections. 

The properties of weak and strong chaos have been observed in many other orbits beyond the four 
examples treated in the paper. After identifying the properties of the q-entropic parameter in typical cases 
of "weakly" and "strongly" chaotic motions, future studies should be made to extend this approach to 
bigger samples of orbits to provide a more quantitative perspective on their distribution in phase and real 
space As explained in [^3], the occurrence of q-Gaussian pdfs allows one to define a new form of partition 
function and thus develop a different type of (nonextensive) thermodynamics to describe dynamical systems 
in these regimes. 

We believe that the above conclusions and results are not limited to barred galaxy models. It would be 
interesting, therefore, to apply the statistical methods described in this paper to different galaxy structures 
to quantify their chaotic regimes in more detail and obtain a more complete picture of their dynamics. 
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